This notebook computes the landcover areas for each point in 'sample_points' based on a 'landcover' vector map and its defined landcover types. The computed areas are exported to a csv file.
This notebook uses GRASS GIS (7.0.4), and must be run inside of a GRASS environment (start the jupyter notebook server from the GRASS command line).
The code in this notebook was modified from examples in "How to write a Python GRASS GIS 7 addon" developed for a FOSS4G Europe 2015 workshop. The code is contained in the python-grass-addon repository on GitHub.com.
landcover – vector map representing land cover
sample_points – vector map containing observer points
radius – radius in map units around each sample point to calculate land cover areas within
use_viewshed – consider only 'visible' land cover, based on a viewshed from each sample point
landcover_area_table_filename — path to export land cover area table as a csv file
n_types — number of land cover types
In [ ]:
landcover = 'landcover'
In [ ]:
landcover_type_field = 'type'
In [ ]:
n_types = 15
In [ ]:
sample_points = 'sample_points_field'
In [ ]:
radius = 500
In [ ]:
use_viewshed = True
get viewshed suffix (for filenames)
In [ ]:
def getViewshedSuffix():
if use_viewshed:
viewshed_suffix = '_viewshed'
else:
viewshed_suffix = ''
return viewshed_suffix
In [ ]:
landcover_area_table_filename = "/home/ubuntu/landcover_areas_{0}m{1}.csv".format(radius, getViewshedSuffix())
In [ ]:
import numpy as np
import pandas
import pyprind
GRASS import statements
In [ ]:
import grass.script as gscript
In [ ]:
from grass.pygrass.vector.geometry import Point
from grass.pygrass.vector import Vector
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.table import DBlinks
connect to attribute table
In [ ]:
def connectToAttributeTable(map):
vector = VectorTopo(map)
vector.open(mode='r')
dblinks = DBlinks(vector.c_mapinfo)
link = dblinks[0]
return link.table()
extract point from vector map
In [ ]:
def extractPoint(input, ID, output):
where = 'ID = {0}'.format(ID)
type = 'point'
gscript.read_command('v.extract',
input=input,
where=where,
output=output,
type=type,
overwrite=True)
create buffer around point
In [ ]:
def bufferPoint(input, output, radius):
gscript.read_command('v.buffer',
input=input,
output=output,
type='point',
distance=radius,
overwrite=True)
convert viewshed from raster to vector map
In [ ]:
def vectorizeViewshed(input, ID, output):
type = 'area'
column = 'visible'
gscript.read_command('r.to.vect',
input=input,
output=output,
type=type,
column=column,
overwrite=True)
overlay a vector map on an underlying vector map using 'and' selection operator
In [ ]:
def overlay(overlay, underlay, output):
operator='and'
gscript.read_command('v.overlay',
ainput=overlay,
binput=underlay,
operator=operator,
output=output,
overwrite=True)
add area column to vector map
In [ ]:
def calculateAreas(map):
#add new area column
gscript.read_command('v.db.addcolumn',
map=map,
columns="area_square_meters DOUBLE PRECISION")
#compute area and insert into area column
gscript.read_command('v.to.db',
map=map,
type='centroid',
option='area',
columns='area_square_meters',
unit='meters')
create table showing total area by landcover type
In [ ]:
def getLandcoverAreaByType(map):
#get area data
table = connectToAttributeTable(map=map)
table.filters.select()
columns = table.columns.names()
cursor = table.execute()
result = np.array(cursor.fetchall())
cursor.close()
data = pandas.DataFrame(result, columns=columns).set_index('cat')
#sum areas by landcover type
data['area_square_meters'] = pandas.to_numeric(data['area_square_meters'])
areas = data[['b_' + landcover_type_field, 'area_square_meters']].groupby(by='b_' + landcover_type_field).sum()
return areas
export vector map to a shapefile
In [ ]:
def exportVectorToShapefile(map, output):
gscript.read_command('v.out.ogr',
input=map,
format='ESRI_Shapefile',
output=output,
flags='e',
overwrite=True)
get info about a vector map
In [ ]:
def getVectorMapInfo(map):
return gscript.read_command('v.info', map=map)
In [ ]:
#view(['1_viewshed'])
# TODO determine how to display vector map
connect to 'sample_points' attribute table
In [ ]:
point_table = connectToAttributeTable(map=sample_points)
point_table.filters.select()
columns = point_table.columns.names()
cursor = point_table.execute()
result = np.array(cursor.fetchall())
cursor.close()
point_data = pandas.DataFrame(result, columns=columns).set_index('cat')
loop through sample points
In [ ]:
with Vector(sample_points, mode='r') as points:
#setup progress bar
progress_bar = pyprind.ProgBar(points.n_lines, bar_char='█', title='Landcover analysis progress', monitor=True, stream=1, width=50)
#iterate through points
for point in points:
#get point ID (SiteID)
ID = point_data['ID'][point.cat-1]
#update progress bar
progress_bar.update(item_id=ID)
#buffer current point
extractPoint(input='sample_points_field', ID=ID, output='tmp_buffer_point')
bufferPoint(input='tmp_buffer_point', output='tmp_point_buffer', radius=radius)
#set buffer as overlay
overlay_input = 'tmp_point_buffer'
#consider only visible landcover if 'use_viewshed' = True
if use_viewshed:
viewshed = 'vect_{0}_viewshed'.format(ID)
visible_viewshed = 'vect_{0}_viewshed_{1}m'.format(ID, radius)
#vectorize viewshed
vectorizeViewshed(input='{0}_viewshed'.format(ID), ID=ID, output=viewshed)
#overlay buffer on viewshed
overlay(overlay=overlay_input,
underlay=viewshed,
output=visible_viewshed)
#set overlay to the visible viewshed
overlay_input = visible_viewshed
overlay_output = 'vect_{0}_landcover_{1}m{2}'.format(ID, radius, getViewshedSuffix())
#overlay landcover
overlay(overlay=overlay_input,
underlay=landcover,
output=overlay_output)
#calculate landcover area
calculateAreas(map=overlay_output)
create table with landcover areas by type for each point
In [ ]:
#create table
index_start = 41
''' set the first index number,
allowing easier insertion into a database table that already contains
area calculations with other parameters
(i.e radius and use_viewshed)'''
columns = ['ID', 'SiteID', 'IncludedArea']
columns = columns + [ str(n) for n in range(1,n_types+1) ]
area_table = pandas.DataFrame(columns=columns)
#set naming variables
included_area = '{0}m{1}'.format(radius, getViewshedSuffix()) #0=radius, 1=viewshed_suffix
map_pattern = 'vect_{0}_landcover_{1}m{2}' #0=ID, 1=radius, 2=viewshed_suffix
#iterate through points
for index, point in point_data.iterrows():
ID = point['ID']
map = map_pattern.format(ID, radius, getViewshedSuffix())
#initiate row
row = {'ID':"{0:.3g}".format(int(index) + index_start),
'SiteID': str(ID),
'IncludedArea': included_area}
#get landcover areas
areas = getLandcoverAreaByType(map)
#iterate through area types
for index, area in areas.iterrows():
#add area to row
row["{0:.3g}".format(int(index))] = area['area_square_meters']
#append row to table
area_table = area_table.append(row, ignore_index=True)
area_table.set_index('ID', inplace=True)
#export table to file
area_table.to_csv(landcover_area_table_filename, header=False)